EE4211 Project¶
Team 1: Quadratic Done by Choi Jae Hyung, Lee Min Young, Nicholas Lui Ming Yang, Tran Duy AnhIn [1]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import datetime
from scipy import stats
from functools import reduce
import sklearn.linear_model as lm
import sklearn
from sklearn.linear_model import LinearRegression as LR
from sklearn.svm import SVR
In [2]:
df = pd.read_csv("project_data.csv")
# df.shape - (1584823, 3)
# df column values:
# localminute - Timestamp
# dataid - MeterID
# meter_value - meter reading
Q1.1¶
How many houses are included in the measurement study?¶
In [3]:
# Since dataid is unique for each meter, we can count the number of unique dataid numbers
df.dataid.value_counts().size
Out[3]:
This gives us 157 for the number of houses included in the study¶
Are there any malfunctioning meters? If so, identify them and the time periods where they were malfunctioning.¶
There are various ways to define a malfunctioning meter. Let us explore some of them¶
Let us first check, if there are meters that report a lower value at a later timestamp. Since meter_value is cumulative consumption, meter_value should not be lower at a later timestamp for the same meter¶
In [4]:
grouped = df.groupby(['dataid'], sort=['localminute'])
def has_decreasing_values(df):
current_value = 0
for index, val in df.iteritems():
if val < current_value:
return True
current_value = val
meters_with_decreasing = (grouped['meter_value']
.agg(has_decreasing_values)
.where(lambda x: x)
.dropna()
.keys())
In [5]:
print(len(meters_with_decreasing))
meters_with_decreasing
Out[5]:
Wow, we have 43 meters that have a decreasing value! Let's zoom in and get the offending sections for each meter. (In other words, find the exact data points which show this descreasing value)¶
In [6]:
## Need to re-sort after filter since filter takes all rows and breaks original sorting
dec_meters = (grouped.filter(lambda x: int(x.name) in meters_with_decreasing)
.groupby(['dataid'], sort=['localminute']))
## Iterate over values to find offending rows for each meter
## WARNING: RUNS VERY SLOWLY. TODO: OPTIMIZE
offending_values = {}
for group_id, rows in dec_meters.groups.items():
offending_values[group_id] = []
current_value = 0
group_rows = dec_meters.get_group(group_id)
group_rows_count = group_rows.shape[0]
for i in range(group_rows_count):
if group_rows.iloc[i]['meter_value'] < current_value:
offending_values[group_id].append([group_rows.iloc[i-1], group_rows.iloc[i]])
current_value = group_rows.iloc[i]['meter_value']
Number of 'broken data' instances for each meter¶
In [7]:
print("Meter ID |", "Number of broken readings")
for k, v in offending_values.items():
print(str(k).ljust(20), len(v))
Just knowing the number of broken readings may not be useful.¶
Let's say that we want to know which meters should be fixed, since faulty meters result in us inaccurately measuring consumption, and possibly losing money.¶
There should be some measure of tolerance in deciding if a meter is broken. In this case, let's check the average decreasing amount, and the ratio of "broken" readings¶
In [8]:
print("Meter ID |", "Number of broken readings |", "Average decrease across broken readings |", "Percentage of broken readings for meter")
for k, v in offending_values.items():
print(str(k).ljust(20),
str(len(v)).ljust(20),
str(reduce(lambda x, y: x + abs(y[1]['meter_value'] - y[0]['meter_value']) / len(v), [0] + v )).ljust(40),
(len(v) / dec_meters.get_group(k).shape[0]) * 100)
In [9]:
print("Meter ID |", "Number of broken readings |", "Average decrease across broken readings |", "Percentage of broken readings for meter")
for k, v in offending_values.items():
measure = str(reduce(lambda x, y: x + abs(y[1]['meter_value'] - y[0]['meter_value']) / len(v), [0] + v )).ljust(40)
if float(measure) > 10:
print(str(k).ljust(20),
str(len(v)).ljust(20),
measure,
(len(v) / dec_meters.get_group(k).shape[0]) * 100)
In [10]:
"""Use this function to validate/check each bad meter given our heuristic"""
def print_bad_meter_readings(meterID):
meter_readings = offending_values[meterID]
print("time apart |".ljust(20), "meter value initial |", "meter value after |", "difference")
for readings_pair in meter_readings:
print(str(pd.to_datetime(readings_pair[1]['localminute']) - pd.to_datetime(readings_pair[0]['localminute'])).ljust(25),
str(readings_pair[0]['meter_value']).ljust(20),
str(readings_pair[1]['meter_value']).ljust(20),
readings_pair[1]['meter_value'] - readings_pair[0]['meter_value'])
print_bad_meter_readings(1556)
Another requirement of the meters is that they push their readings to be saved if the meter values differ by at least 2 cubic foot¶
Let us verify that every pair of readings for a meter differ by at least 2 cubic foot. To not double count, we will only check for readings where the differing value is 0 >= x > 2 so that we don't get the same error readings where the readings decrease¶
Since the value is smaller, we will not focus on the difference in values, but on the percentage of total readings for the meter only¶
In [11]:
grouped2 = df.groupby(['dataid'], sort=['localminute'])
def has_stagnant_values(df):
current_value = 0
for index, val in df.iteritems():
if index == 0:
current_value = val
continue
if val < current_value + 2 and val >= current_value:
return True
current_value = val
meters_with_stagnant = (grouped['meter_value']
.agg(has_stagnant_values)
.where(lambda x: x)
.dropna()
.keys())
print("Number of meters with stagnant values: ", len(meters_with_stagnant))
The following segment should have calculated the percentage of stagnant values for each meter, but it takes too long (tried for 5 minutes and gave up). At least we know number of meters that reported stagnant values?¶
In [12]:
# ## Need to re-sort after filter since filter takes all rows and breaks original sorting
# stagnant_meters = (grouped2.filter(lambda x: int(x.name) in meters_with_stagnant)
# .groupby(['dataid'], sort=['localminute']))
# ## Iterate over values to count offending occurences. Not counting value
# ## WARNING: RUNS VERY SLOWLY. TODO: OPTIMIZE
# offending_values2 = {}
# for group_id, rows in stagnant_meters.groups.items():
# offending_values2[group_id] = 0
# group_rows = stagnant_meters.get_group(group_id)
# group_rows_count = group_rows.shape[0]
# # Set current_value so we do not trigger first reading
# current_value = group_rows.iloc[0]['meter_value'] - 5
# for i in range(group_rows_count):
# if group_rows.iloc[i]['meter_value'] < current_value + 2 and group_rows.iloc[i]['meter_value'] >= current_value:
# offending_values2[group_id] += 1
# current_value = group_rows.iloc[i]['meter_value']
# print("Meter ID |", "Number of broken readings |", "Percentage of broken readings for meter")
# for k, v in offending_values2.items():
# print(str(k).ljust(20),
# str(v).ljust(20),
# (v / stagnant_meters.get_group(k).shape[0]) * 100)
In [13]:
stagnant_meters = (grouped2.filter(lambda x: int(x.name) in meters_with_stagnant)
.groupby(['dataid'], sort=['localminute']))
print("Total number of readings to iterate, which is why it is slow: ",
reduce(lambda x, y: x + len(y), [0] + list(stagnant_meters.groups.values())))
We now have two different criterion for deciding what constitues a broken meter, and the readings that helped us determine that. Since the question asks for "time periods where they were malfunctioning.", let's demonstrate that we can consolidate broken readings into time periods¶
In [14]:
# let's use the meter with the most broken readings (based off decreasing values)
most_broken_values_meter = grouped.get_group(5403)
In [15]:
# We are using a stricter requirement now. Let's consider that, if the subsequent meter reading is less
# than an increment of 2 (meters are only supposed to report after at least an increment of 2)
# Broken_criteria is a helper function (comparator function) that takes in two (in-order) readings and outputs a
# boolean of whether the later reading is "broken"
broken_criteria = lambda x, y: y['meter_value'] < x['meter_value'] + 2
In [16]:
broken_readings = 0
num_readings = most_broken_values_meter.shape[0]
for i in range(1, num_readings):
if broken_criteria(most_broken_values_meter.iloc[i-1], most_broken_values_meter.iloc[i]):
broken_readings += 1
print("Number of broken readings for meter 5403, based on stricter requirements: ", broken_readings)
In [17]:
# Let's now create a function that can aggregate a meter's readings into "broken" time periods
# Return value of the function will be as follow:
"""
2D array.
First dimension is the time periods.
Second dimension is the consecutive broken readings that make up the broken time period
To find the actual time data, take the ['localminute'] attribute from the first and last reading per period
time_periods = [
[reading_1, reading_2, reading_3],
[reading_1, reading_2]
]
"""
def get_broken_time_periods(broken_criteria, meter_readings):
num_readings = meter_readings.shape[0]
time_periods = []
temp_period = []
for i in range(1, num_readings):
if broken_criteria(meter_readings.iloc[i-1], meter_readings.iloc[i]):
temp_period.append(meter_readings.iloc[i])
else:
if temp_period:
time_periods.append(temp_period)
temp_period = []
if temp_period:
time_periods.append(temp_period)
return time_periods
In [18]:
broken_time_periods = get_broken_time_periods(broken_criteria, most_broken_values_meter)
print("Number of broken time periods: ", len(broken_time_periods))
Objective:¶
Obtain hourly data from a list of given data, and plot the obtained data
Assumption (How we aggregate the data into hourly readings)¶
- The data is always cumulative: increasing in value as time goes.
The lowest value in an hour should be the closest value to an exact hour (hh:00:00) and a good estimate of that hour's data.
Therefore, the first data point in an hour is used to represent the hourly data. (The assumption is that we are interested in approximating the consumption amount at the start of the hour)
How to handle missing data?¶
As the data is cumulative gas consumption by household, when there is no data, it can be assumed as there is insignificant gas consumption between that period.
Hence, the missing data will be given same value as most recent hourly data obtained.
Bad data readings¶
This plot of hourly readings does not account for "bad" data readings (consecutive reading with a lower value, which should not happen for cumulative consumption data)
The likelihood of seeing an hourly data point that is bad is Average P(bad readings an hour) (460), since we have max 4 readings per minute
This reduces the code complexity to get a plot, and the usefulness of this approach depends on the usage of this plot
If we are simply looking for a general consumption level visualisation, then this plot should suffice. However, if we are looking for plots with a decrease in value to identify faulty meters, then this plot is insufficient. For that, we cannot hold our original assumption (that readings are cumulative and strictly increasing), and should take the MIN(hour_values) to represent that hour
Code¶
The algorithm of hourly data collection consists of the functions below.¶
- In every new hour, collect the lowest value as the hourly data.
- Look out for missing hour. When missing hour is found, previuos hour data is given for that hour.
Hourly Data Collection¶
Create DateHour column, that represents "Date + Hour" of the data.¶
In [6]:
#change to datetime format
df["localminute"] = pd.to_datetime(df["localminute"])
#get a dataframe of data by hour:
df_byhour = pd.DataFrame(df.groupby([df["localminute"].dt.date.rename("Date"),
df["localminute"].dt.hour.rename("Hour"),
df["dataid"]]).min().reset_index())
# To preview df_byhour sorted by meter and time
# df_byhour = df_byhour.sort_values=(by=['dataid', 'Date', 'Hour'])
#add a column called "DateHour" in the dataframe:
df_byhour['DateHour'] = pd.to_datetime(df_byhour['Date'].apply(str)+' '+df_byhour['Hour'].apply(str)+':00:00')
Choose the meter_ID, starting time and ending time of the hourly data collection.¶
In [7]:
#Format: ID No.
ID = 739
#ranges from 06/10/2015 00:00:00 to 05/11/2015 23:00:00
start_date = datetime.date(2015, 10, 5)
end_date = datetime.date(2015, 11, 6)
filteredDF = df_byhour[(df_byhour["Date"] > start_date) & (df_byhour["Date"] < end_date)]
Collect hourly data.¶
In [8]:
#from 06/10/2015 -> 06/11/2015 there are 31 days -> 744 data points
prevHr = pd.Timestamp('2015-10-06 00:00:00')
endHr = pd.Timestamp('2015-11-05 23:00:00')
#generate a temporary dataframe through get_group function with each id
dfTemp = filteredDF.groupby(['dataid']).get_group(ID)
#get the list of hours and meter values
hourList = dfTemp['DateHour'].tolist()
valList = dfTemp['meter_value'].tolist()
#re-initialize the prev val as the first value of the val list
prevVal = valList[0]
#re-initialize the list of values as a list with 1 item
newList = [valList[0]]
#loop through each hour in the list
for index, hr in enumerate(hourList):
#if hour is more than an hour bigger than previous
while (hr - prevHr).seconds >= 3600:
newList.append(prevVal)
#increment every hour
prevHr += pd.Timedelta(seconds=3600)
#update the prev value everytime one index is passed
prevVal = valList[index]
#some data does not reach until the end date, thats where this comes in
while prevHr < endHr:
newList.append(prevVal)
#increment every hour
prevHr += pd.Timedelta(seconds=3600)
Plot the hourly reading.¶
In [12]:
x = pd.date_range(datetime.date(2015, 10, 6), datetime.date(2015, 11, 6), freq = 'H').tolist()
#get rid of the last hour which is 00:00:00 06/11/2016
x.pop()
y = newList
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
ax.plot(x, y, 'ok', ms=1)
plt.xlabel ("Time")
plt.ylabel ("Meter value")
ax.set_title('Hourly Data')
plt.show()
In [13]:
#get a list of meter id
idList = df['dataid'].unique()
#there are 183 days => 4392 hours => 4392 readings in each list
#generate a dictionary
hourlyDataDict = {}
#define start hour and end hour
staHr = pd.Timestamp('2015-10-01 05:00:00')
endHr = pd.Timestamp('2016-04-01 04:00:00')
#loop through each id
for id1 in idList:
#re-innitialize the prev Hour as start hour
prevHr = staHr
#generate a temporary dataframe through get_group function with each id
dfTemp = df_byhour.groupby(['dataid']).get_group(id1)
#get the list of hours and meter values
hourList = dfTemp['DateHour'].tolist()
valList = dfTemp['meter_value'].tolist()
#re-initialize the prev pal as the first value of the val list
prevVal = valList[0]
#re-initialize the list of values as a list with 1 item
newList = [valList[0]]
#loop through each hour in the list
for index, hr in enumerate(hourList):
#if hour is more than an hour bigger than previous
while (hr - prevHr).seconds >= 3600:
newList.append(prevVal)
#increment every hour
prevHr += pd.Timedelta(seconds=3600)
#update the prev value everytime one index is passed
prevVal = valList[index]
#some data does not reach until the end date, thats where this comes in
while prevHr < endHr:
newList.append(prevVal)
#increment every hour
prevHr += pd.Timedelta(seconds=3600)
#after all that create a new entry in the Dictionary with the key as the id and value is the list of meter readings
hourlyDataDict[id1] = newList
Creating a Dictionary of Correlation Coefficient.¶
In [14]:
corDict = {}
#loop through id List to get the first id
for id1 in idList:
#create the empty value for the first key
corDict[id1] = {}
#loop through the id list to get the second id
for id2 in idList:
# generate the coefficient through numpy.corrcoef (can be changed later)
coef = np.corrcoef(hourlyDataDict[id1], hourlyDataDict[id2])[0, 1]
# assign the value to the dict with the appropriate key
corDict[id1][id2] = coef
#print("coefficient between " + str(id1) + " and " + str(id2) + " is: " + str(coef))
Now, let's find the top 5 houses with highest correlation value for each house¶
In [15]:
from collections import Counter
idList.sort()
top5dict = {}
for id in idList:
#### This is if you just want the houseID for each (first column being the ID)
# print(sorted(corDict[id], key =corDict[id].get, reverse=True)[:6])
#### This is if you want to get the houseID + the value of correlation coefficient
c = Counter(corDict[id])
mc = c.most_common(6)
del mc[0]
top5dict.update({id : mc})
top5dict
Out[15]:
Q1.4 Student Proposal¶
Focus: Find out which meter needs maintenance¶
- Draw histogram to see the number of faulty meters with various magnitude of decreasing value. (Prof. Mehul's recommendation to visualize)
On top of magnitude, no. of times the value decrease can be second histogram to observe.
- Set our own standard from the visualization and calculate number of meters in need of maintenance
Follow-up idea:¶
1) Print a summary to inform the company & consumer about the meter's status. "We recommend maintenance because __."
2) Expected amount of payment for the maintenance. "Expected price: $__. Please dial 1800-900-XXXX"
Q2.1 In this part, you will asked to build a model to forecast the hourly readings in the future (next hour). Can you explain why you may want to forecast the gas consumption in the future?¶
Predicting gas consumption in the future is necessary because of many reasons:
- Gas is a natural resource and it is limited, predicting gas consumption will shed some light onto how much more we need to produce and when it is necessary to cut down on the consumption for fear of exhausting this natural resource.
- It allows the gas companies to purchase in bulk in the anticipation of the consumption. This could potentially reduce purchase price, cost of transferring gas, and inventory fee.
- Since forecast is based on the current consumption pattern, if the real data deviates from forecast by certain percentage, then it would be easier to investigate what has caused the deviation from forecast. i.e. Environment campaign, unexpected cold wind drift causing people to use more heater, etc.
Who would find this information valuable?¶
Just to list down a few: customers, government, environmental activists, gas company, maintenance company, gas production countries.
What can you do if you have a good forecasting model?¶
With a good forecasting model, there could be various actions that can be done: Sell the data / models to
- Gas companies
- Maintenance companies
- Government or environmental activitists
This might be also able to predict
- When should next servicing of the system be done before it becomes faulty.
- Show trends for amount of gas consumed in a region or consumption in coming months or even years. Then, gas companies can prepare in advance for bulk orders for cheaper cost of production. Or could reduce inventory issues by preparing only the suitable amount of gas.
- Possible effect of consumption of natural gases on global warming in the future.
Q2.2 Build linear regression to forecast the hourly readings in the future.¶
(Regarding Q2.2 and Q2.3 future reading estimation beyond the 6 months data, they are shown at the end of Q2.3)
In [16]:
hourArray = features_train = np.asarray(range(4392)).reshape(-1,1)
clfDictLR = {}
for key in hourlyDataDict:
#each key is a linear regression model,
clfDictLR[key] = LR()
#features_train = the number of hour from 5:00:00 01/10/2015,
features_train = hourArray
#print(features_train),
#labels_train = values,
labels_train = np.asarray(hourlyDataDict[key])
#fit the model accordingly,
clfDictLR[key].fit(features_train, labels_train)
hourArray = features_train = np.asarray(range(4392)).reshape(-1,1)
In [17]:
staHr = pd.Timestamp('2015-10-01 05:00:00')
dtStar = datetime.datetime(2015, 10, 1) + datetime.timedelta(hours=5)
dtEnd = datetime.datetime(2016, 4, 1) + datetime.timedelta(hours=4)
#get a list of dateTime from the start to the end\n
dateTimeList = pd.date_range(dtStar, dtEnd, freq = 'H').tolist()
hourList = [[((x - staHr).seconds + (x - staHr).days*86400)/ 3600] for x in dateTimeList]
def newPredLR(meter_id):
return clfDictLR[meter_id].predict(hourList)
i) Time Series Plot¶
In [18]:
for key in hourlyDataDict:
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
plt.xlabel('hours')
plt.ylabel('values')
x = dateTimeList
y = np.asarray(hourlyDataDict[key])
ax.plot(x, y,color='g')
ax.plot(x, newPredLR(key) ,color='k')
plt.title("Linear Regression Time Series Plot, Meter ID = " + str(key))
plt.show()
ii) Scatter Plot¶
In [19]:
for key in hourlyDataDict:
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
plt.xlabel('hours')
plt.ylabel('values')
x = dateTimeList
y = np.asarray(hourlyDataDict[key])
ax.scatter(x, y,color='g',s=0.05)
ax.plot(x, newPredLR(key) ,color='k')
plt.title("Linear Regression Scatter Plot, Meter ID = " + str(key))
plt.show()
Q2.3: Build Support Vector regression to forecast the hourly readings in the future.¶
Achieving more accurate prediction model through averaging the hourly values between all 157 meters.¶
As shown in Q2.2, there are multiple meters that are having malfunctioning. Amongst the malfunctioning meters, some jump in meter readings exceptionally. This spike increase may cause the prediction to be inaccurate.
Assumption:¶
There is a similar pattern in household gas consumptions. Hence taking the average value of all the gas consumption of households every hour gives general hourly consumption thereby reducing the estimation error made by huge jump in values of some meters.
Therefore, by obtaining average hourly consumption value, this general consumption is used to predict the future gas consumption.¶
1) Modify the hourlyDataDict to hourlyConsumptionDataDict¶
In [20]:
hourlyConsumptionDataDict = {}
for key in hourlyDataDict:
# all values in hourly data dict is subtraced by the first value
hourlyConsumptionDataDict[key] = [x - hourlyDataDict[key][0] for x in hourlyDataDict[key]]
# calculate the average across all meter for each hour
avgDataList = []
for i in range(4392):
total = 0
for key in hourlyConsumptionDataDict:
total += hourlyConsumptionDataDict[key][i]
#print('total is: ' + str(total))
avg = total/157
avgDataList.append(avg)
2) Build Linear Regression and Support Vector Regression classifier and fit the data¶
In [21]:
hourArray = features_train = np.asarray(range(4392)).reshape(-1,1)
features_train = hourArray
labels_train = np.asarray(avgDataList)
clfLR = LR()
clfLR.fit(features_train, labels_train)
clfSVR = SVR(kernel='linear')
clfSVR.fit(features_train, labels_train)
Out[21]:
3) Generate a list of datetime values for predictions¶
In [22]:
staHr = pd.Timestamp('2015-10-01 05:00:00')
dtStar = datetime.datetime(2015, 10, 1) + datetime.timedelta(hours=5)
dtEnd = datetime.datetime(2016, 4, 1) + datetime.timedelta(hours=4)
#get a list of dateTime from the start to the end\n
dateTimeList = pd.date_range(dtStar, dtEnd, freq = 'H').tolist()
hourList = [[((x - staHr).seconds + (x - staHr).days*86400)/ 3600] for x in dateTimeList]
In [23]:
staHr = pd.Timestamp('2015-10-01 05:00:00')
dtStar_future = datetime.datetime(2015, 10, 1) + datetime.timedelta(hours=5)
dtEnd_future = datetime.datetime(2016, 10, 1) + datetime.timedelta(hours=4)
#get a list of dateTime from the start to the end\n
dateTimeList_future = pd.date_range(dtStar_future, dtEnd_future, freq = 'H').tolist()
hourList_future = [[((x - staHr).seconds + (x - staHr).days*86400)/ 3600] for x in dateTimeList_future]
4) Plot the new Linear Regression and Support Vector Regression on the same plot for comparison¶
i) Time Series Plot¶
In [24]:
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
plt.xlabel('hours')
plt.ylabel('values')
x = dateTimeList
y = np.asarray(avgDataList)
ax.plot(x, y,color='g')
ax.plot(x, clfLR.predict(hourList) ,color='r')
ax.plot(x, clfSVR.predict(hourList) ,color='b')
plt.title("Linear Regression and Support Vector Regression, with Times Series Plot")
plt.show()
ii) Scatter Plot¶
In [25]:
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
plt.xlabel('hours')
plt.ylabel('values')
x = dateTimeList
y = np.asarray(avgDataList)
ax.scatter(x, y,color='g',s=0.05)
ax.plot(x, clfLR.predict(hourList) ,color='r')
ax.plot(x, clfSVR.predict(hourList) ,color='b')
plt.title("Linear Regression and Support Vector Regression, with Scatter Plot")
plt.show()
Effect of averaging the data¶
By averaging the data, the consumption generally increases, looking more like a proper cumulative gas consumption data.
Is Linear Regression & Support Vector Regression the same?¶
Linear Regression & Support Vector Regression looks almost the same!
Let's find out by predicting the future value, beyond the 6 months data given.
Future Data Prediction¶
Time series plot is shown for future data prediction.
In [26]:
fig, ax = plt.subplots(1, 1, figsize=(20, 6))
plt.xlabel('hours')
plt.ylabel('values')
x = dateTimeList
y = np.asarray(avgDataList)
ax.plot(x, y,color='g')
ax.plot(dateTimeList_future, clfLR.predict(hourList_future) ,color='r')
ax.plot(dateTimeList_future, clfSVR.predict(hourList_future) ,color='b')
plt.title("Predicting Future Data with Linear Regression and Support Vector Regression, with Times Series Plot")
plt.show()
From the above future data prediction plot, we can see that the linear regression line and support vector regression line is actually not the same!¶
Q3.1 Student Proposal¶
As per section 1.4:
Focus: Find out which meter needs maintenance
Draw histogram to see the number of faulty meters with various magnitude of decreasing value. (Prof. Mehul's recommendation to visualize)
On top of magnitude, no. of times the value decrease can be second histogram to observe.
Set our own standard from the visualization and calculate number of meters in need of maintenance
Follow-up idea:
1) Print a summary to inform the company & consumer about the meter's status. "We recommend maintenance because __."
2) Expected amount of payment for the maintenance. "Expected price: $__. Please dial 1800-900-XXXX"
In [180]:
# offending_values is the variable storing all "broken meter" values
# That is consecutive readings with decreasing consumption
# Let's format this data for plotting a histogram of
# magnitude of decreasing values, to determine if there is a pattern to the broken meters
# (is there a systematic decrease)?
decreasing_histogram_values = []
for readings in offending_values.values():
for reading in readings:
decreasing_histogram_values.append({
"meter_id": reading[0].dataid,
"decrease_magnitude": reading[0].meter_value - reading[1].meter_value
})
decreasing_values_df = pd.DataFrame(decreasing_histogram_values)
Histogram of decreasing values for all meters¶
In [182]:
decreasing_values_hist = decreasing_values_df.hist(bins=100, column=['decrease_magnitude'])
In [184]:
broken_reading_count_df = decreasing_values_df.groupby(['meter_id']).agg(lambda x: x.shape[0]).rename({'decrease_magnitude': 'broken_count'}, axis='columns')
# Add in meters with no broken readings
appending_values = []
offending_keys = offending_values.keys()
for key in grouped.groups.keys():
if key not in offending_keys:
appending_values.append({
"meter_id": key,
"broken_count": 0
})
broken_reading_count_all_df = broken_reading_count_df.append(pd.DataFrame(appending_values).set_index('meter_id'))
broken_reading_count_all_hist = broken_reading_count_all_df.hist(bins=5)
In [185]:
broken_reading_count_hist = broken_reading_count_df.hist(bins=50)
In [186]:
broken_reading_count_filtered_hist = (broken_reading_count_df
.where(broken_reading_count_df['broken_count'] > 5)
.dropna()
.hist(bins=50))
With this, it is clear that most meters have only a few (<5) broken readings, or do not actually have any broken readings.¶
However, this may not be very useful, as some meters may not have reported as many readings¶
(less values reported results in less "chances" to report a broken reading)¶
Let's fix this by drawing a histogram of frequency of broken readings for all meters¶
In [187]:
broken_freq_values = []
for key in grouped.groups.keys():
broken_values = 0
broken_freq = 0
total_values = grouped.get_group(key).shape[0]
if key in offending_values:
broken_values = len(offending_values[key])
broken_freq = broken_values / total_values * 100
broken_freq_values.append({
"meter_id": key,
"broken_readings": broken_values,
"total_readings": total_values,
"broken_frequency": broken_freq
})
broken_freq_df = pd.DataFrame(broken_freq_values).set_index('meter_id')
broken_freq_hist = broken_freq_df.hist(bins=10, column="broken_frequency")
In [188]:
broken_freq_filtered_hist = (broken_freq_df
.where(broken_freq_df['broken_frequency'] > 0)
.hist(bins = 10, column="broken_frequency"))
In [189]:
broken_readings_cumsum = broken_reading_count_all_df.sort_values(['broken_count']).cumsum()
# change index beceause meter_id is no longer in-order, and then plot cumulative sum
broken_readings_cumsum.reset_index().plot(y="broken_count", title="cumulative sum of broken readings")
Out[189]:
In [190]:
### Additional stats to support hypothesis
print ("70% of broken readings: ", 0.7 * broken_readings_cumsum.max().broken_count)
print ("Number of meters required to fill 70% of readings: ",
broken_readings_cumsum
.where(broken_readings_cumsum.broken_count
> 0.7 * broken_readings_cumsum.max().broken_count).dropna().shape[0])
How does this apply to real life considerations?¶
There are several factors that affect the ability and desire of a company to repair a meter¶
Here are some assumptions that we will make, to simplify our modelling process¶
- There are resources dedicated towards repairing meters
- If unused, they are simply wasted
- As such, the problem becomes one of resource allocation. Where can we direct our resources to maximise their effectiveness?
We have earlier explored this as 3 separate measures
1) Absolute number of broken readings
2) Frequency of broken readings
3) Magnitude (or severity) of broken readings
A scoring system can be implemented, to give each meter a "broken" score across a certain time window. This will then be used to prioritize resources
What if the company does not have dedicated resources for repairing meters?¶
- The assumption now is that there is a cost incurred for fixing meters, since there is no dedicated resource to do so
- Thus, the problem is now one of calculating the loss incurred by a meter, from the meter readings
- Same as above, we are relying on the same 3 measures
These factors are related because they can be aggregated to calculate the losses incurred from each meter owing to erroneous readings.
What is missing from this data set is the actual gas consumption. With that, we can be sure of the actual extent of misreporting incurred be erroneous readings
Is there anything else we can realise from looking at our data? Let's look at the ten most broken meters¶
In [191]:
# 8156, 5403
ten_most_broken_meter_ids = broken_readings_cumsum.tail(10).index.values
for meter_id in ten_most_broken_meter_ids:
grouped.get_group(meter_id).plot(x="localminute", y="meter_value", title=str(meter_id)+" consumption")
Breakthrough!¶
It is clear that most of these broken readings occured during the same timeframe for these meters.¶
Rather than marking these meters as broken and wasting time fixing them, we can instead seek to understand what caused the broken readings.
- The actual company will need to investigate likely causes during the error timeframe to determine this
Instead of simply marking decreasing readings to determine whether meters are broken, we can additionally check if the decreasing readings are instead a result of spikes, and whether they normalize after some time